function idat = fsb_interpolate_parallel_scans(idat,sandbox)

% FSB : Interpolate image data for artefact periods
%
% EXAMPLE:
% idat = fsb_interpolate_parallel_scans(idat,sandbox)
%
% INPUT:
% idat: 4-D image data
% sandbox: sandbox experiment struct
%
% OUTPUT:
% idat: 4-D image data with interpolated artefact periods
%
% CALLED BY:
% FSB.m
%
% NOTES:
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if correct dataset size
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

d = size(idat);
if size(d)<4;
    disp('No 4D dataset, aborting');
    return
end

if sandbox.fullscan == 0;
    disp('Concatenated dataset, aborting');
    return
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run function with sandbox input
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

explicit = 0;

for scannum = 1:max(sandbox.intrial(:,4));

    scanindex = find(sandbox.intrial(:,4) ==scannum);
    trialnum = length(unique(sandbox.intrial(scanindex,3)))-1;

    trialind = find(sandbox.intrial(scanindex,2)>0 & sandbox.intrial(scanindex,1)~=0);
    dummy1 = zeros(length(scanindex),1);
    dummy1(trialind) = 1;
    dummytest = bwlabeln(dummy1);

    for volnum = 1:trialnum;
        trialind2 = find(dummytest==volnum);
        testi(volnum) = length(trialind2);
    end

    otto = mean(testi);

    idat_seg = single(idat(:,:,:,scanindex));

    %~~~~~~~~~~~~~~~~~~~~~~~~
    % Create indices of trial beginning and end
    %~~~~~~~~~~~~~~~~~~~~~~~~
    trialindex1 = find(sandbox.intrial(scanindex,2)<3 & sandbox.intrial(scanindex,1)~=0);
    trialindex2 = find(sandbox.intrial(scanindex,2)>otto-2 & sandbox.intrial(scanindex,1)~=0);
    trialindexdummy1 = zeros(length(scanindex),1);
    trialindexdummy1(trialindex1) = 1;
    trialindexdummy2 = zeros(length(scanindex),1);
    trialindexdummy2(trialindex2) = 1;
    test1 = bwlabeln(trialindexdummy1);
    test2 = bwlabeln(trialindexdummy2);
    g = waitbar(0);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Go through trials and figure out what needs to be interpolated
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    for volnum = 1:trialnum;

        detstring = ['Calculation artefact interpolation for trial ' num2str(volnum) ' of ' num2str(trialnum)];
        disp(detstring)
        waitbar(volnum/trialnum,g,detstring);
        volind1{volnum} = find(test1==volnum);
        startvol = idat_seg(:,:,:,volind1{volnum});
        startvolt{volnum} = mean(startvol,4);
        volind2 = find(test2==volnum);
        endvol = idat_seg(:,:,:,volind2);
        endvolt{volnum} = mean(endvol,4);

    end
    close(g);
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Now interpolate it
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    g = waitbar(0);
    if trialnum>1
        for volnum = 2:trialnum;

            detstring = ['Interpolating artefacts of trial ' num2str(volnum) ' of ' num2str(trialnum)];
            disp(detstring)
            waitbar(volnum/trialnum,g,detstring);
            voldiff = startvolt{volnum} - endvolt{volnum-1};
            volind2 = find(test2==volnum-1);

            try
                artlength = volind1{volnum}(1)-max(volind2);
            catch
                disp(volind2);
            end

            artsec = repmat(endvolt{volnum-1},[1 1 1 artlength-1]);
            voldiff2 = voldiff/artlength;
            voldiff2 = repmat(voldiff2,[1 1 1 artlength-1]);
            artgrad2 = 1:artlength-1;
            artgrad2 = repmat(artgrad2,[size(voldiff2,1) 1 size(voldiff2,2) size(voldiff2,3)]);
            artgrad2 = permute(artgrad2,[1 3 4 2]);
            artgrad3 = artgrad2.*voldiff2;
            artgrad3 = artgrad3+artsec;
            try
                idat_seg(:,:,:,max(volind2)+1:volind1{volnum}(1)-1) = artgrad3;
            catch
                dispi
            end

        end

    else
        volnum =1;
    end

    try
        endind = find(test2==volnum);
        endseg = max(endind)+1:size(idat_seg,4);
        idat_seg(:,:,:,endseg) = repmat(endvolt{volnum},[1 1 1 length(endseg)]);
    catch
        disp('Issues with end segment')
        endind = [];
    end
    startseg = 1:min(volind1{1})-1;

    idat_seg(:,:,:,startseg) = repmat(startvolt{1},[1 1 1 length(startseg)]);

    disp(['...done with scan ' num2str(scannum)]);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Calculate some values for display purposes if desired
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if explicit ==1;
        trialchange(scannum) = sigchange2(scannum)/trialnum;
        volchange(scannum) = sigchange2(scannum)/length(scanindex);
        voltrial = max(sandbox.intrial(scanindex,2));
        trialpart = (voltrial*trialnum)/length(scanindex);
        trialchange(scannum) = trialchange(scannum)*trialpart;
    end;

    idat(:,:,:,scanindex) = idat_seg;

    close(g);
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create plots
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if explicit ==1;
    figure(501);
    set(gcf,'Name',' Plot of mean volume intensity before (blue) and after (red) Detrending');

    plot (avgint); hold on;
    plot (avgint2,'r');

    figure(502);
    set(gcf,'Name',' Plot of mean signal intensity change');
    subplot(1,2,1);
    bar(sigchange2);
    axis tight
    xlabel('Run');
    ylabel('% Signal change per run');
    subplot(1,2,2);
    bar(trialchange);
    axis tight
    xlabel('Trial');
    ylabel('% Signal change per trial');

    figure(503);
    set(gcf,'Name',' Plot of mean signal intensity change');
    plot(avgint3);
end

end
